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1.0  SUMMARY 


In  this  study,  the  tasks  of  correlating  observations  and  detecting  maneuvers  of  space 
objects  are  addressed  through  application  of  minimum- fuel  distance  metrics,  which  provide  a 
natural  tool  with  which  to  associate  space  object  observation  data.  A  numerical  trajectory 
optimization  software  is  developed  to  solve  deterministic  minimum-fuel  boundary  value 
problems,  specifically  for  impulsive  and  low-thrust  trajectories.  A  framework,  referred  to  as  the 
Data  Association  using  Minimum-Fuel  Maneuver  Metric  (DAMMM)  algorithm,  is  implemented 
to  utilize  the  trajectory  optimization  software  and  control  distance  metrics  for  maneuver 
hypothesis  testing,  incorporating  existing  theory  accounting  for  boundary  condition  uncertainty. 
A  more  detailed  explanation  of  this  algorithm  and  the  components  therein  can  be  found  in 
Appendix  A.  A  simulated,  realistic  scenario  is  constructed  to  test  the  algorithm,  and  the  results 
are  compared  to  real  data  from  a  similar  scenario  gathered  from  Wide  Area  Augmentation 
System  (WAAS)  data.  Using  the  simulated  and  operational  data  results,  the  algorithm  is  tested  to 
evaluate  its  limitations  and  tuned  for  application  to  operational  use.  The  control  distance  metric 
data  association  technique  is  compared  to  other  data  association  techniques  (e.g.  Euclidean  and 
Mahalonobis  distances),  revealing  the  merits  of  the  control  distance  metric  in  orbital  mechanics 
applications. 

2.0  INTRODUCTION 

The  problem  of  properly  correlating  on-orbit  observations  and  detecting  object 
maneuvers  is  a  challenging  and  persistent  endeavor.  There  are  currently  at  least  17,000  trackable 
on-orbit  objects,  1,000  of  which  are  active  [1,  2],  and  these  numbers  are  expected  to  grow 
significantly  due  to  improved  tracking  capabilities,  new  launches,  and  continued  debris 
generation  [3].  Predicting  conjunction  events  is  a  difficult  task  [4];  however,  recent  events 
highlight  the  mutual  interest  that  national  and  private  operators  share  for  accurate  object 
correlation  and  maneuver  detection  capability  [5]. 

The  objective  of  this  study  is  to  validate  a  novel  method  for  correlating  space  object 
tracks  with  known  Resident  Space  Objects  (RSOs)  and  to  characterize  and  reconstruct  RSO 
maneuvers  in  support  of  the  Space  Situational  Awareness  (SSA)  Mission  Area.  The  problem  of 
correlating  Uncorrelated  Tracks  (UCTs)  over  large  time  periods  is  particularly  difficult  when 
objects  maneuver  during  observation  gaps.  Even  relatively  small  station-keeping  maneuvers  at 
Geostationary  Earth  Orbit  (GEO)  can  result  in  position  discrepancies  of  many  kilometers  after  an 
observation  gap.  UCT  correlation  is  further  confounded  by  state  estimate  uncertainties.  Since 
both  the  initial  and  final  UCTs  are  best  estimates,  with  associated  uncertainty  distributions,  the 
question  of  correlation  becomes  difficult  to  answer  in  operational  settings,  particularly  in 
densely-populated  regions  of  the  space  environment. 

Given  a  propagated  best  estimate  of  the  state  and  its  associated  distribution,  correlating 
UCTs  tests  the  hypothesis  that  a  newly  observed  state  (with  its  associated  distribution)  could 
possibly  be  a  previously  observed  object,  and  if  not,  what  the  discrepancy  is.  There  are  many 
distance-  or  pseudo-distance  metrics  that  may  be  used  to  measure  the  distance  between  two  state 
distributions  (e.g.  Mahalonobis  distance).  Problematically,  none  of  these  distance  metrics 
directly  quantify  the  level  of  propulsive  effort  required  to  effect  the  observed  state  change.  This 
difference  is  critical,  as  very  small  fuel  expenditures  at  specific  points  in  an  orbit  can  produce 
outsized  state  discrepancies. 
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This  study  uses  a  minimum-fuel  control  distance  metric  and  hypothesis  testing  to 
approach  data  association,  maneuver  detection  and  characterization  problems.  Because  on-board 
fuel  remains  a  scarce  commodity  for  operational  spacecraft,  operators  are  likely  to  execute 
optimal  or  near-optimal  maneuvers.  Under  the  assumptions  of  optimal  control,  multiple 
deterministic  observations  can  be  related  by  computing  the  control  effort  required  for  a 
spacecraft  to  meet  those  boundary  conditions.  This  approach  necessitates  the  reconstruction  of  a 
minimum- fuel  trajectory  consistent  with  a-priori  information  and  new  observations,  supporting 
SSA  maneuver  characterization  needs.  Optimal  connecting  trajectories  are  computed  in  a 
MATLAB  minimum  quadratic  control  trajectory  optimization  implementation  of  trajectory 
optimization  software.  The  resulting  optimal  trajectory  is  used  as  an  input  in  the  proposed 
hypothesis  testing  framework,  incorporating  anomaly  detection  results  from  the  previous  papers 
[6,  7,  8].  The  combined  framework  is  tested  using  simulated  and  real-world  operational  data  to 
validate  the  method  and  evaluate  any  limitations.  Finally,  the  results  from  these  simulation  cases 
are  compared  to  data  association  approaches  using  other  distance  metrics. 


2.1  Statement  of  Work 

The  original  Statement  of  Work  (SoW)  for  this  effort  is  enumerated  below.  Each  SoW  element 
is  explicitly  identified  in  the  appropriate  sections  of  this  report.  Minor  modifications  to  the  SoW 
agreed  upon  by  the  Program  Manager  and  Principal  Investigator  (PI)  are  identified  and  described 
in  footnotes  where  appropriate. 

1.  Framework  Implementation.  Incorporate  existing  theory  accounting  for  boundary 
condition  uncertainty  and  anomaly  detection  hypothesis  testing.  In  particular,  the 
prototype  code  used  to  generate  results  in  previous  literature  will  be  rewritten  to 
accommodate  more  general  test  cases  and  data  inputs.' 

2.  Compare  Results.  Compare  the  control  distance  metric  data  association  techniques 
against  other  data  association  approaches  (e.g.,  Mahalanobis  distance,  orbit  element 
gating)  using  simulated  data^.  After  the  proposed  approach  has  been  implemented  and 
internally  tested,  the  results  will  be  compared  with  previously  proposed  and  currently 
used  approaches  such  as  Mahalonobis  distance  and  orbit  element  gating.  As  the  project 
progresses,  different  data  association  methods  will  be  compared. 

3.  Extend  Theory^.  The  theoretical  approach  to  reducing  variational  cost  conservatism  will 
be  pursued  using  analytical  methods.  After  a  satisfactory  theory  accounting  for  boundary 
condition  variations  has  been  developed,  it  will  be  compared  with  both  the  more 


'  For  programmatic  purposes,  the  development,  implementation,  and  inclusion  of  trajectory 
optimization  code  was  also  added  to  the  effort. 

^  During  the  performance  period  of  this  effort  the  empirical  WAAS  data  became  available  and  a 
similar  analysis  with  this  data  is  included  in  this  report. 

The  Galaxy  15  trajectory  information  in  the  WAAS  data  presumes  an  electric  thruster 
propulsion  system;  trajectory  fuel  cost  function  modifications  were  unnecessary  for  this  effort. 
An  extension  of  the  theory  to  allow  for  Gaussian  Sum  representations  of  the  boundary  constraint 
uncertainties  is  included  instead  in  Appendix  B. 
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conservative  quadratic  cost  upper  bound  and  the  unseented  transform  approaeh  recently 
developed  in  the  literature. 

The  remainder  of  this  report  addresses  these  SOW  items  as  follows:  SoW  item  1,  ‘Framework 
Implementation,’  is  discussed  in  seetion  2.2,  DAMMM  Overview,  SoW  item  2,  ‘Compare 
Results,’  is  addressed  in  section  3.0  and  4.0,  and  SoW  item  3,  ‘Extend  Theory,’  is  diseussed  in 
Appendix  B. 

2.2  DAMMM  Overview 

To  address  this  researeh  goal,  the  algorithm  framework  (and  MATLAB  software  block 
diagram)  is  defined  in  Figure  1  and  deseribed  below.  The  variables  of  Figure  1  are  as  follows: 

^  is  the  mean  state  of  the  observation  at  t^,  g  is  mean  state  of  the  observation  at  tg, 
PxA  is  the  covariance  of  the  observation  at  Px,b^^  eovariance  of  the  observation  tg, 
X  (t)  is  the  optimal  state  trajeetory,  u  (t)  is  the  optimal  eontrol  trajectory,  p  (t)  is  the  optimal 
costate  trajectory,  is  the  quadratie  eontrol  eost  of  the  optimal  trajectory,  f(J)  is  the  quadratie 
control  cost  distribution,  and  P^j  is  the  probability  of  anomaly.  This  software  framework  and 
associated  eode  documentation  in  Appendix  A  constitute  the  work  satisfying  the  ‘Framework 
Implementation’  SoW  item  in  section  2.1.  The  inputs  to  the  algorithm  are  two  UCTs,  containing 
the  estimated  states  of  the  trajeetory  boundary  eonditions  along  with  covarianee  matrices 
deseribing  the  probability  distribution  of  these  boundary  states.  The  first  portion  of  the  algorithm 
simply  parses  the  input  traeks  to  extraet  the  nominal  boundary  states  at  the  initial  point  (labeled 
as  point  A  in  Figure  1)  and  final  point  (B)  and  stores  the  eovariance  data  for  use  later. 


f(J) 

Figure  1.  DAMMM  Algorithm  Functional  Flow  Diagram 


The  nominal  states  are  used  to  generate  an  optimal  trajeetory  that  conneets  the  boundary 
points  through  the  trajectory  optimization  algorithm.  This  algorithm  solves  the  deterministic  two- 
point  boundary  value  problem  using  MATLAB’s  eonstrained  minimization  function,  fmincon. 
The  eost  function  being  minimized,  shown  in  Equation  1,  is  the  quadratic  control  cost  (J),  which 
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is  a  function  of  the  control  thrust  accelerations  (u)  and  is  ideal  for  variable  speeifie  impulse 
engines  (VSI): 


J  =  ‘/2lu^udt  (1) 

The  initial  and  final  points  are  used  as  boundary  eonditions  and  the  simulation  time  is  diseretized 
into  a  user-defined  number  of  intermediate  steps.  The  deeision  variable  for  this  minimization  is  a 
staeked  vector  of  the  thrust  aeceleration  vectors  at  each  discrete  time  instant.  Keplerian  two-body 
dynamies  are  enforced  between  steps  of  the  trajectory  as  equality  eonstraints  to  ensure  the 
generated  trajectory  dynamics  are  aecurate.  Beeause  the  partial  derivatives  of  the  dynamics  with 
respect  to  the  deeision  variables  (thrust  aceelerations)  are  well  known,  the  gradient  of  the 
constraint  is  supplied  to  the  optimization  algorithm,  greatly  improving  performance  of  the 
implementation.  The  output  of  this  step  is  a  nominal  optimal  state  and  control  trajectory 
connecting  the  UCTs.  The  generated  optimal  trajeetory  is  validated  using  the  nonlinear  dynamics 
to  numerically  integrate  the  proposed  control  vector  and  quantify  the  error  between  the 
integrated  final  condition  and  the  speeified  final  UCT  boundary  condition.  Adjusting  tunable 
convergence  parameters  of  the  optimization  algorithm,  sueh  as  absolute  and  relative  toleranees, 
can  reduce  this  error.  For  this  study,  the  relative  and  absolute  convergence  criteria  are  both  set  to 
10'*^  to  decrease  final  state  error.  The  well-behaved  nature  of  the  dynamies  allows  solutions  to 
even  this  low  toleranee  on  the  order  of  tens  of  steps,  so  the  effect  on  computational  tractability  is 
minimal.  For  instance,  over  100  test  runs  of  the  nominal  simulated  trajeetory  (defined  in  Section 
3.1),  the  trajeetory  optimization  routine  ran  5  iterations  of  the  optimizer  before  converging,  with 
an  average  runtime  of  0.314  seeonds.  The  eomputer  used  for  this  test  is  an  i5-2500K  clocked  at 
4.2  GHz  with  16  GB  random  aecess  memory  (RAM)  running  64-bit  MATLAB  2013a. 

The  optimal  states  x  (t)  and  controls  u  (t)  for  Ia  <  t  <  te  are  then  used  to  estimate  a  history 
of  optimal  eontrol  eostates  p  (t).  This  step  is  required  because  algorithm  that  computes  the 
control  cost  distribution  requires  costates  instead  of  thrust  accelerations.  The  costates  are 
constructed  iteratively  by  guessing  a  eostate  trajeetory  based  on  the  optimal  eontrol  trajeetory 
and  then  modifying  the  costates  through  iterative  least-squares  method  until  the  eontrol  cost  is 
within  a  user-defined  eonvergence  toleranee  of  the  optimal  trajeetory.  As  an  additional  cheek, 
the  new  optimal  costate  trajeetory  is  compared  to  the  original  optimal  trajeetory  to  ensure  the 
endpoint  error  has  not  ehanged  signifieantly. 

Using  the  optimal  state  and  costate  trajeetories,  as  well  as  the  covariance  information 
from  the  input  traeks,  control  probability  distributions  are  generated  using  the  techniques 
outlined  in  a  previous  paper  [8].  Finally,  the  nominal  optimal  control  cost  and  control 
distribution  are  used  to  test  the  hypothesis  that  a  maneuver  actually  occurred  and  eomplete  the 
maneuver  deteetion  task. 

3.0  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 

As  described  in  seetion  Seetion  2.2,  eaeh  UCT  pair  eombination  examined  in  this  study  is 
considered  a  Two-point  Boundary  Value  Problem  (TPBVP).  The  trajeetory  optimization  is 
performed  to  generate  a  nominal  optimal  trajeetory  conneeting  the  two  boundary  eonditions 
(UCTs).  Presently,  Keplerian  two-body  dynamics  are  assumed  and  other  perturbations  are 
ignored  (e.g.,  lunar,  solar,  non-spherieal  perturbations).  For  each  case  tested,  a  nominal  UCT 
boundary  condition  uncertainty  (covariance  matrix)  is  parameterized  using  a  scalar  multiplier  to 
examine  the  effects  of  uncertainty  on  anomaly  detection.  Similarly,  to  examine  the  effects  of 
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time  separation  between  UTCs,  each  scenario  simulates  various  time  intervals  between  UCTs  of 
up  to  48  hours.  The  resulting  control  cost  distributions  due  to  state  uncertainty  at  the  boundary 
conditions  are  approximated  using  the  analytic  first  and  second  moments  as  normal  distributions 
as  specified  in  Holzinger  et  al  [8],  The  deterministic  control  cost  is  then  compared  to  the  control 
cost  distributions  to  rigorously  test  the  maneuver  /  anomaly  hypothesis. 

Using  these  boundary  conditions,  a  deterministic  maneuver  connecting  UCT  A  and  UCT 
B  is  computed.  After  generating  the  nominal  optimal  control  trajectory  and  the  corresponding 
costate  trajectory  (x  (t)  and  p  (t)),  uncertainties  in  position  and  velocity  are  used  to  generate  an 
approximate  probability  distribution  of  the  control  cost.  The  position  and  velocity  uncertainty 
values  used  in  this  simulation  are  shown  in  Table  1.  The  mean  of  the  total  cost  distribution  (np)  is 
equal  to  the  optimal  nominal  control  cost  (/)  plus  the  minimum  detectable  cost  for  a  given 
uncertainty  condition,  where  the  minimum  detectable  cost  is  calculated  from  the  state  uncertainty 
boundary  conditions.  The  process  for  calculating  the  mean  (/Tp)  and  standard  deviation  (Cp)  from 
the  optimal  trajectory  and  covariance  is  based  on  Equations  2  and  3  respectively.  The  details  of 
the  derivation  of  these  equations  are  in  the  previous  paper  [8]  with  Tr  being  the  trace  operator. 

Mp=r  +  Tr[nPJ  (2) 

(7p  =  oi^PzO)  +  2Tr[ftP2ftP2]  (3) 

Here,  Pz  is  a  12-by-12  matrix  with  the  initial  and  final  covariances  on  the  upper-left  and  lower- 
right  6x6  block  diagonals,  Q  is  the  term  in  the  performance  function  corresponding  to  the 
quadratic  in  the  perturbation  term  (Equation  15  in  Reference  8),  and  w  is  the  term  in  the 
performance  function  corresponding  to  the  cross-term  between  the  nominal  control  distance  and 
perturbation  (Equation  14  in  Reference  8).  Furthermore,  the  minimum  detectable  quadratic  cost 
is  used  to  determine  if  state  uncertainty  accounts  for  the  state  change,  and  is  defined  as  the  mean 
of  the  cost  distribution  generated  from  state  uncertainty,  shown  in  Equation  4  (using  terms  from 
Equation  2): 


Jaet  =  Tr[np,]  (4) 

In  order  to  better  model  the  cost  distribution,  an  exact  distribution  is  also  constructed  using  a 
Monte-Carlo  approach  using  a  random  vector  SZ  G  N(0,  P^).  The  quadratic  cost  due  to  this  state 
uncertainty  can  then  be  calculated  as  shown  in  Equation  5  [8]: 

June  =  +  SZ'^SISZ  (5) 

Ten-thousand  costs  were  sampled  and  binned  into  1000  bins  to  construct  a  distribution  of  costs. 
The  minimum  detectable  cost  for  this  distribution  is  then  found  by  computing  the  normalized 
Probability  Density  Function  (PDF)  and  Cumulative  Distribution  Function  (CDF).  The  cost 
where  the  CDF  for  this  control  distribution  equals  0.5  is  defined  as  the  minimum  detectable 
control  cost.  While  both  Gaussian  approximate  and  Monte-Carlo  distributions  are  generated,  the 
resulting  distributions  are  nearly  identical,  so  the  remainder  of  the  report  uses  the  (more 
accurate)  Monte-Carlo  distribution. 

A  trajectory  with  a  nominal  control  cost  J  greater  than  the  minimum  detectable  cost  J^et 
has  a  higher  probability  (greater  than  0.5,  asymptotically  approaching  1  as  J  increases)  of  being 
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the  result  of  a  controlled  maneuver  or  unmodeled  perturbation  and  not  just  the  result  of 
propagated  boundary  condition  uncertainty. 

The  probability  that  a  hypothesis  is  true  (p^)  can  be  computed  from  the  distribution  of 
costs  (//•(/*))  as  shown  in  Equation  6  [7]: 

Pa  =  L'l  /y-aOdy*  (6) 

A  sample  result  from  this  method  (using  the  synthetic  scenario  defined  in  Section  3.1)  is  shown 
in  Figure  3.  The  figure  shows  the  PDF  and  CDF  for  the  quadratic  control  cost  as  a  function  of 
control  cost.  The  maximum  value  of  the  normalized  PDF  corresponds  to  the  minimum  detectable 
cost  as  defined  in  Equation  4,  and  the  vertical  dashed  line  represents  the  nominal  optimal  control 
cost  for  that  maneuver.  On  the  CDF  plot,  it  can  be  seen  that  the  intersection  of  the  CDF  plot  and 
nominal  quadratic  cost  yields  a  probability  of  anomaly  of  approximately  0.75,  or  75%.  Indeed, 
the  CDF  is  the  key  to  computing  the  probability  that  a  maneuver  has  occurred  for  a  given 
optimal  trajectory  and  state  uncertainty.  This  probability  of  a  maneuver  occurring  is  also  referred 
to  as  the  probability  of  anomaly,  denoting  the  possibility  that  anomalous  behavior  outside  the 
assumed  dynamics  (e.g.  a  maneuver  or  unmodeled  perturbations)  are  the  cause  for  the  state 
change. 


Figure  2.  Normalized  PDF  and  CDF  vs  Difference  From  Minimum  Detectable  Control  Cost 


Since  the  control  cost  and  distribution  varies  for  each  time  horizon  and  UCT 
uncertainties,  it  is  more  intuitive  to  define  and  use  a  Signal-to-Noise  Ratio  (SNR)  as  the  ratio  of 
the  nominal  optimal  control  cost  to  the  minimum  detectable  cost,  as  shown  in  Equation  7: 


T* 

SNR  = 

Jdet 


J* 

Tr[nPz] 


(7) 
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If  the  nominal  optimal  control  cost  (the  signal)  is  equal  to  the  minimum  detectable  quadratic  cost 
for  a  given  state  uncertainty  (the  noise),  then  the  SNR  is  unity.  When  SNR  =  1 ,  the  probability  of 
anomaly  is  equal  to  0.5,  or  50%.  That  is  to  say,  there  is  an  equal  probability  that  the  observed 
state  change  is  caused  by  an  anomaly  (such  as  a  control  maneuver)  or  caused  by  uncertainty  in 
the  initial  and  final  states.  As  the  nominal  cost  increases  with  respect  to  the  minimum  detectable 
cost,  the  SNR  increases  and  the  probability  that  the  state  change  is  caused  by  an  anomaly 
increases.  As  SNR  increases,  the  nominal  cost  becomes  substantially  larger  than  the  cost 
uncertainty,  so  the  probability  that  state  uncertainty  could  cause  the  state  change  is  zero  and  the 
probability  of  anomaly  approaches  1.0,  or  100%. 

Given  an  optimal  trajectory  and  uncertainty  in  the  boundary  conditions,  one  can  compute 
the  likelihood  that  the  observed  state  change  (or  lack  thereof)  is  due  to  quiescent  propagation  or 
active  maneuvering,  thus  addressing  the  maneuver  detection  goal.  This  final  step  in  the 
DAMMM  procedure  is  referred  to  as  testing  the  anomaly  hypothesis,  as  shown  in  Figure  1. 
Additionally,  if  it  is  determined  that  boundary  condition  uncertainty  cannot  account  for  the  state 
change,  and  a  maneuver  or  unmodeled  perturbation  must  have  occurred,  the  minimum  magnitude 
and  time  history  of  the  potential  maneuver  can  be  computed,  characterizing  any  such  maneuver. 
Furthermore,  if  the  maneuver  is  determined  to  be  impossible  given  the  assumed  propulsive 
capability  of  the  RSO,  this  method  can  be  used  to  associate  (or  disassociate,  in  this  case)  UCTs 
based  on  the  control  authority  required. 

The  first  part  of  SoW  item  2,  ‘Compare  Results,’  is  to  validate  the  software 
implementation.  To  do  so,  extreme  test  cases  are  used  to  identify  the  limits  and  performance  of 
this  approach  in  non-ideal  situations.  Primarily,  the  process  sensitivity  to  boundary  condition 
uncertainty  and  simulation  duration  must  be  understood  to  better  utilize  and  interpret  the  results 
of  the  algorithm.  Increasing  state  uncertainty  decreases  the  SNR  (minimum  detectable  cost 
increases),  and  it  is  important  to  know  what  SNRs  yield  reliable  results.  Additionally,  simulation 
duration  (the  time  between  initial  and  final  boundary  conditions)  affects  the  performance  of  the 
trajectory  optimizer,  and  the  problem  of  correlating  UCTs  over  large  time  periods  is  particularly 
difficult  when  objects  maneuver  during  gaps  in  observations. 

To  validate  the  software  implementation,  a  sensitivity  study  is  conducted  to  generate 
confidence  curves  of  probability  of  anomaly  based  on  a)  maneuver  SNR  and  b)  UCT  time 
separation.  The  boundary  condition  uncertainty  is  initialized  at  1  meter  in  position  and  0.01 
meters-per-second  in  velocity,  as  shown  in  Equation  8: 

Pa  =  Pb  —  diag(a^,  a^,  a^,  (O.Ola)^,  (O.Ola)^,  (O.Ola)^)  ,  (8) 

where  a  is  the  scalar  that  is  varied  between  scenarios  to  change  the  uncertainty  (and  modulate  the 
effective  SNR).  The  scalar  was  allowed  to  vary  logarithmically  from  1  to  10^,  determined  by  trial 
and  error  to  cover  a  wide  range  of  SNR  values  (from  near-zero  to  above  6).  This  corresponds  to 
position  and  velocity  uncertainties  of  100  km  /  1  km/s  at  maximum,  and  10m  /  O.lm/s  at 
minimum.  The  simulation  duration  between  UCTs  is  varied  as  well,  keeping  the  ascending  node 
passage  at  the  midpoint  of  the  simulation  duration.  The  initial  and  final  times  are  varied  from  5 
minutes  before  and  after  the  ascending  node  respectively  (10  minute  simulation  duration)  to  24 
hours  before  and  after  the  ascending  node  respectively  (48  hour  simulation  duration).  The  results 
of  this  study,  as  well  as  the  nominal  1-hour  duration  case,  are  included  in  Section  4.1. 
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3.1  Scenario  Descriptions 


To  validate  the  DAMMM  framework,  two  scenarios  are  considered.  First,  a  synthetic 
data  (simulated)  scenario  is  constructed  with  a  geostationary  spacecraft  executing  a  one-degree 
inclination-change  maneuver  at  the  ascending  node.  This  particular  scenario  was  chosen  to 
replicate  the  available  WAAS  data  for  the  Galaxy  15  geostationary  spacecraft.  The  second 
scenario  uses  the  WAAS  data,  which  enables  examination  of  Galaxy  15  N-S  (inclination)  electric 
thruster  station-keeping  maneuvers.  This  process  illustrates  a  possible  operational  use-case  of  the 
software  package  and  identifies  how  the  proposed  DAMMM  framework  for  maneuver  detection, 
characterization,  and  UCT  correlation  can  work.  Each  scenario  is  described  individually  in  the 
following  subsections. 

3.1.1  Synthetic  Data  Scenario. 

The  synthetic  (simulated)  scenario  is  constructed  to  emulate  an  inclination-change  station¬ 
keeping  maneuver,  such  as  the  operational  case  provided  in  the  WAAS  data.  The  boundary 
conditions  for  this  simulation  can  be  seen  in  Table  1.  A  spacecraft  is  initialized  in  a  circular  orbit 
at  GEO  with  an  initial  inclination  of  one-degree.  The  final  state  has  an  identical  semi-major  axis, 
eccentricity,  and  right  ascension  of  the  ascending  node,  now  has  zero  inclination.  Note  that,  since 
this  orbit  is  circular,  the  argument  of  periapsis  is  undefined.  Notice  also  that  the  true  anomaly 
changes  by  roughly  15  degrees  during  the  simulation.  At  GEO,  this  approximately  corresponds 
to  a  one-hour  simulation  duration.  The  initial  condition  (UCT  A)  occurs  30-minutes  before  the 
spacecraft  crosses  the  ascending  node  (where  the  maneuver  is  executed),  and  the  final  condition 
(UCT  B)  occurs  30-minutes  after  the  ascending  node  crossing.  A  time-step  of  one  minute 
between  each  discrete  step  of  the  simulation  is  selected,  so  the  simulation  is  divided  into  61 
equal  steps.  Figure  3  shows  the  inclination  approaching  the  target  value  over  the  course  of  the 
simulation,  as  well  as  the  commanded  control  (thrust  acceleration),  in  meters-per-second. 


Table  1.  Simulated  Scenario  Boundary  Conditions 


Parameter 

Initial  Condition 

Final  Condition 

Units 

(UCT  A) 

(UCT  B) 

Semi-major  Axis 

42164 

42164 

kilometers 

Inclination 

1 

0 

degrees 

Eccentricity 

0 

0 

- 

Right  Ascension  of 
Ascending  Node 

0 

0 

degrees 

Argument  of 
Periapsis 

n/a 

n/a 

degrees 

True  Longitude 

352.4795 

7.5205 

degrees 

Position  Uncertainty 

la 

la 

meters 

Velocity 

Uncertainty 

O.Ola 

O.Ola 

meters-per- 

second 
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Figure  3.  Sample  Trajectory  and  Control  for  Plane  Change  Maneuver 


The  synthetic  data  scenario  is  used  in  two  analyses  in  this  study.  The  first  is  the  specific 
test  case  detailed  above,  with  varying  state  uncertainty,  testing  the  outcome  of  the  anomaly 
hypothesis  process.  The  second  is  a  sensitivity  study,  meant  to  extend  the  algorithm  to  longer 
duration  simulations  and  more  extreme  SNR  situations,  to  test  the  algorithm  capabilities. 


3.1.2  Real-World  Data 

To  complement  the  simulated  scenario,  the  algorithm  is  tested  using  real  operational  data 
supplied  by  Captain.  Kulumani  and  Dr.  Baldwin  at  the  Air  Force  Research  Laboratory  (AFRL) 
Space  Vehicles  Directorate.  The  availability  of  this  real-world  data  drove  the  construction  of  the 
simulated  scenario  in  Section  3.1.1.  The  data,  taken  from  observations  of  the  Galaxy  15 
geostationary  satellite,  spans  the  month  of  February  2014,  and  includes  Earth-centered  Earth- 
fixed  (ECEF)  position  and  velocity,  as  well  as  radial,  in-track,  and  cross-track  (RIC)  acceleration 
(as  seen  by  a  rotating  Hill  frame  attached  to  the  spacecraft).  Inspection  of  the  acceleration  data, 
shown  in  Figure  4,  reveals  two  candidates  for  the  North-South  (inclination)  station-keeping 
maneuver:  one  during  February  9*,  and  another  during  February  23'^'*.  Both  candidates  show  a 
large  anomalous  cross-track  acceleration  (seen  in  the  bottom  plot  of  Figure  4),  which  is 
indicative  of  a  plane  change  maneuver  and  therefore  a  station-keeping  maneuver. 
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Figure  4.  RIC  Acceleration  Data  of  Galaxy  15  from  WAAS  Data 


The  data  is  analyzed  in  a  similar  manner  to  the  synthetic  data  sensitivity  study:  by 
varying  boundary  condition  uncertainty  and  time,  a  sensitivity  study  is  conducted  to  construct 
probability  of  anomaly  confidence  curves  with  respect  to  SNR  and  simulation  duration.  From  the 
two  candidate  maneuvers,  the  initial  boundary  condition  is  selected  as  the  February  8*'^ 
maneuver,  just  before  the  large  cross-track  thrust.  This  initial  boundary  condition  is  held  fixed, 
and  the  final  state  and  simulation  duration  are  changed  by  choosing  different  final  boundary 
condition  states.  Therefore,  the  study  begins  with  a  simulation  duration  of  1  hour  and  increasing 
to  a  duration  of  48  hours.  The  results  of  this  sensitivity  study  are  included  in  Section  4.2. 

4.0  RESULTS  AND  DISCUSSION 

The  results  from  this  effort  are  broken  down  into  numerous  sections,  addressing  various 
portions  of  the  proposed  work.  The  results  from  the  simulated  scenario  are  discussed  to  validate 
the  minimum-fuel  distance  metric  approach  to  maneuver  detection.  Similarly,  the  real-world  data 
is  used  to  test  the  application  of  the  software  suite  to  a  plausible  operational  situation  (part  1  of 
SoW  item  ‘Compare  Results’).  The  minimum-fuel  distance  metric  approach  is  also  compared  to 
other  common  distance  metrics  to  assess  its  usefulness,  particularly  in  spacecraft  in  maneuver 
detection  (part  2  of  SoW  item  ‘Compare  Results’). 
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4.1  Framework  Validation 


The  simulated  one-hour,  one-degree  GEO  inclination  correction  maneuver,  as  defined  in 
Table  1  and  Section  3.1,  is  used  to  validate  the  DAMMM  algorithm,  ensuring  functionality  of  the 
algorithm  as  shown  in  Figure  1.  The  resulting  nominal  simulated  maneuver  trajectory  is  shown 
in  orbit-element  space  in  Figure  2,  showing  the  1 -degree  inclination  correction.  The  quadratic 
cost  for  this  maneuver  is  shown  in  Table  2.  Using  the  control  history  obtained  from  the  trajectory 
optimizer,  the  actual  velocity  change  (denoted  as  AV)  can  be  computed  as  the  integral  of  the 
thrust  accelerations  over  time.  As  a  sanity  check,  this  velocity  change  is  compared  to  other 
bounding  costs  for  validation.  For  instance,  any  non-impulsive  maneuver  is  bounded  from  below 
by  the  impulsive  equivalent  of  that  maneuver.  The  velocity  change  required  for  an  impulsive  1- 
degree  inclination  change  at  GEO  is  shown  in  Table  2,  and  optimal  trajectory  is  bounded  from 
below  by  the  impulsive  cost  as  expected.  Additionally,  an  upper-bound  cost  can  be  applied  based 
on  the  quadratic  cost  from  the  maneuver  [9].  The  value  of  the  upper  bound  for  this  scenario  is 
also  shown  in  Table  2.  The  trajectory  cost  is  within  the  acceptable  range  for  both  of  these  checks. 


Table  2.  Simulated  Maneuver  Control  Cost  and  Comparisons 


Parameter 

Value 

Units 

Quadratic  Control  Cost 

0.4022 

m^/s"' 

Impulsive  AV 

53.6624 

m/s 

Actual  AV 

53.8151 

m/s 

Upper-bound  AV 

53.8153 

m/s 

Using  the  nominal  optimal  trajectory  for  this  scenario,  uncertainty  is  introduced  in  the 
initial  and  final  UCTs.  The  initial  and  final  UCTs  are  each  initialized  with  1-meter  position 
uncertainty  and  0.01-meter-per-second  velocity  uncertainty.  Notice  that  the  velocity  uncertainty 
is  much  smaller  than  the  required  velocity  change  from  Table  2,  so  intuition  says  that  this 
condition  is  highly  observable  and  the  probability  of  anomaly  is  1 .  Indeed,  this  is  the  case  as  the 
nominal  SNR  for  this  case  is  very  high.  The  uncertainties  were  then  scaled  together  up  to  lO"^ 
times  the  initial  amount  in  order  to  increase  the  position  and  velocity  noise,  for  a  maximum 
position  and  velocity  uncertainty  of  100  kilometers  and  1  kilometer-per-second,  respectively. 
Figure  5  shows  the  result  of  the  variation  of  the  probability  of  anomaly  as  a  function  of  SNR.  As 
desired,  the  probability  of  anomaly  at  an  SNR  of  1  is  exactly  0.5,  so  there  is  an  equal  likelihood 
that  the  observed  state  change  is  caused  by  a  control  maneuver  or  uncertainty.  As  SNR  increases 
above  1 ,  the  likelihood  that  a  maneuver  occurred  increases,  and  vice  versa  below  an  SNR  of  1 .  It 
can  be  noted  that  above  an  SNR  of  6  the  probability  of  anomaly  is  always  1,  so  the  plots  have 
been  truncated  to  only  show  the  interesting  lower  ranges  of  SNR.  Given  the  results  from  this 
simulated  scenario,  if  the  SNR  is  above  2,  the  algorithm  would  say  that  there  is  a  high  likelihood 
that  either  the  object  seen  has  maneuvered  or  the  tracks  are  not  correlated  and  actually  belong  to 
separate  objects.  In  order  to  make  this  association  distinction,  prior  information  on  the  object  and 
the  propulsive  requirements  from  the  optimal  trajectory  should  be  considered.  This  nominal 
simulated  scenario  has  shown  that  the  DAMMM  algorithm  is  capable  of  ingesting  UCT  data 
(states  and  distributions)  and  utilizing  knowledge  of  orbital  mechanics  to  provide  data  to  aid  in 
the  association  of  the  tracks. 
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Figure  5.  Probability  of  Anomaly  vs.  SNR,  Simulated  Scenario,  1-hour  Duration 


The  algorithm  has  been  shown  to  function  well  in  this  specific  scenario,  but  in  order  to 
better  characterize  the  performance  of  the  algorithm,  a  sensitivity  study  is  conducted.  The 
variable  parameters  are  the  state  uncertainty  (as  before)  and  the  time  between  observations, 
otherwise  known  as  the  simulation  duration,  which  was  fixed  in  the  previous  analysis.  Similar  to 
the  previous  analysis,  the  state  uncertainty  was  initialized  at  1 -meter  position  uncertainty,  0.01- 
meter-per-second  velocity  uncertainty,  and  was  scaled  up  to  10^  times  to  vary  SNR.  The 
uncertainty  had  to  be  varied  so  drastically  because  the  minimum  detectable  quadratic  cost  varies 
widely  with  changes  in  simulation  duration.  Simulation  durations  varied  between  10  minutes  and 
48  hours,  varying  evenly  on  either  side  of  the  ascending  node  crossing. 

4.2  Synthetic  Data  Parametric  Performance 

Figure  6  is  a  contour  plot  containing  the  results  from  this  investigation;  SNR  and  duration 
are  varied  on  the  x-  and  y-axes  respectively,  and  the  contours  show  lines  of  constant  probability 
of  anomaly.  The  contours  shown  represent  different  increasing  probabilities  of  anomaly  that  the 
maneuver  would  be  observed  given  the  observation  duration  and  uncertainties.  Once  again,  the 
0.5,  or  50%,  confidence  level  is  a  straight  line  at  an  SNR  of  1,  as  expected  due  to  the  nature  of 
the  anomaly  detection  method.  The  results  show  a  relatively  low  impact  of  simulation  duration 
on  the  detectability  of  the  maneuver  since  the  confidence  contours  are  mostly  straight  lines.  The 
major  exception  to  this  occurs  at  simulation  durations  of  12  and  36  hours,  which  correspond  to 
the  initial  and  final  conditions  with  the  highest  out-of-equatorial-plane  state  difference  between 
the  observations.  At  these  points  in  the  orbit,  smaller  boundary  condition  uncertainty  is  required 
to  be  able  to  conclude  that  a  maneuver  has  indeed  occurred.  However,  for  even  75%  confidence 
in  maneuver  detection,  the  SNR  required  is  never  greater  than  2. 
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Figure  6.  Anomaly  Detection  Probability  Confidence-level  Contours  for  Uncertainty  and 
Duration  Sensitivity  Study,  Synthetic  Scenario 


This  result  emphasizes  an  advantage  of  the  control  metric  approach,  which  is  that  it 
performs  consistently  with  respect  to  observation  duration.  Some  previous  knowledge  of  the 
orbit  is  still  useful  in  determining  the  validity  of  the  maneuver  detection  probability  results.  For 
instance,  when  attempting  to  detect  inclination  change  maneuvers,  observations  spaced  at  half- 
orbit  periods  are  not  ideal  because  this  is  could  yield  a  large  state  difference  and  require  lower 
observation  uncertainty  to  detect  maneuvers  with  high  likelihood.  Another,  more  obvious,  result 
is  that  larger  maneuvers,  which  yield  a  higher  SNR  for  the  same  boundary  condition  uncertainty, 
are  easier  to  detect.  This  is  particularly  relevant  for  maneuvers  by  spacecraft  with  variable 
specific-impulse  engines,  where  the  state  changes  are  small  and  spread  out  over  a  long  time. 
Conversely,  a  large  state-change  by  a  nearly-impulsive  maneuver  could  still  be  detected  even  in 
a  short  time  between  observations,  if  the  control  required  for  the  maneuver  significantly  exceeds 
the  state  change  in  the  homogenous  orbit  based  on  boundary  condition  uncertainty. 

4.3  WAAS  Data  Parametric  Performance 

To  further  investigate  the  effectiveness  of  this  algorithm,  and  to  apply  it  to  a  situation 
representative  of  operational  situations,  real  data  was  obtained  from  AFRL  and  used  in  the 
algorithm  as  described  in  Section  3.2.  The  data  obtained  was  from  roughly  one  month  of 
observations  (during  February  2014)  of  the  Galaxy  15  spacecraft  in  geostationary  orbit.  This 
spacecraft  features  an  IHI  BT-4  propulsion  system,  which  performs  as  a  variable  specific- 
impulse  engine,  so  the  results  should  be  in  line  with  section  4.1.  The  data  was  examined  to  pick 
out  point  in  time  where  a  station-keeping  maneuver  likely  occurred,  determined  by  the 

13 

Approved  for  public  release;  distribution  is  unlimited. 


acceleration  data  shown  in  Figure  4.  As  discussed  in  Section  3.2,  the  initial  boundary  condition 
was  chosen  as  a  state  just  before  the  large  cross-track  acceleration  observation  on  February  8*. 
Beginning  with  this  point,  the  final  condition  was  selected  by  increasing  the  duration  between 
observations,  up  to  48-hour  simulation  duration.  The  state  uncertainty  was  varied  similar  to  the 
simulated  case,  beginning  with  1 -meter  position  uncertainty  and  0.0 1-meter-per- second  velocity 
uncertainty  and  scaling  up  to  a  factor  of  10^  times  to  vary  SNR. 

The  results  of  this  sensitivity  study  are  shown  in  Figure  7.  The  confidence-level  contours 
show  similarities  with  Figure  6  from  the  synthetic  scenario.  The  0.5  confidence  level  is  once 
again  a  straight  line  at  a  SNR  of  1,  as  prescribed  by  the  method.  In  this  case,  though,  the 
probability  of  detection  is  even  more  consistent.  This  is  likely  because  of  the  differences  between 
the  modeled  dynamics  and  real  system.  The  simulation  assumes  Keplerian  two-body  dynamics, 
whereas  in  reality  many  accelerations  are  acting  on  the  RSO,  including  n-body  gravitational 
effects,  non-spherical  gravitational  perturbations,  and  solar  radiation  pressure.  As  seen  in  the 
acceleration  data  in  Figure  4,  the  spacecraft  is  thrusting  continuously,  likely  rejecting  some  of  the 
perturbations  to  maintain  its  orbit  slot.  Due  to  this  constant  maneuvering  and  unmodeled 
dynamics,  there  is  less  correlation  between  states  with  half-orbit-period  observation  time 
differences  than  in  the  simulation.  Therefore  the  nominal  control  cost  required  to  associate  UCTs 
is  more  consistent  throughout  the  orbit. 

The  WAAS  data  was  also  sampled  with  a  lower  duration  between  observations,  with  an 
average  20-minute  duration  increase  between  the  sampled  WAAS  data  points.  The  agreement 
between  the  simulated  and  operational  data  cases  reinforces  the  applicability  of  this  algorithm  to 
real  data.  Even  given  the  simplifications  in  the  algorithm,  notably  the  Keplerian  two-body 
dynamics  assumptions  of  the  trajectory  optimization  step,  the  results  from  using  real  data  set 
closely  resemble  the  simulated  example.  This  lends  confidence  to  the  application  of  this 
algorithm  to  operational  data. 


Figure  7.  Anomaly  Detection  Probability  Confidence-level  Contours  for  Uncertainty  and 

Duration  Sensitivity  Study,  WAAS  Data 


14 

Approved  for  public  release;  distribution  is  unlimited. 


4.4  Comparisons  to  Other  Distance  Metrics 

The  proof  of  control  performance  as  a  metric  is  shown  in  the  previous  work  by  proving  the 
properties  of  positivity,  strict  positivity,  symmetry,  and  triangle  inequality  [8],  In  order  to  assess 
the  effectiveness  of  the  control  cost  as  a  distance  metric,  it  is  important  to  consider  comparisons 
to  other  common  distance  metrics.  Primarily,  it  is  important  to  note  the  differences  between 
control  cost  and  both  geometric  and  multivariate  distribution  distance  metrics,  and  the  advantage 
that  a  control  cost  metric  lends  in  orbital  mechanics  applications. 


Figure  8.  Geometric  Distance  Metrics  Applied  to  Maneuvered  Object  Correlation 


A  pictorial  example  of  geometric  distance  metrics  is  included  in  Figure  8.  The  most 
common  geometric  distance  metric  for  evaluating  the  difference  between  two  orbital  state 
vectors  is  Euclidean  distance.  This  considers  only  the  position  difference  between  two  different 
states  without  consideration  for  uncertainty  in  the  boundary  conditions.  The  Euclidean 
distance  is  calculated  by  taking  the  2-norm  of  the  difference  in  position,  expressed  in  Equation  9: 

=  \  \[l  0]{Xg  —  Xj^(t^^^\\2  ,  (9) 

where  xb  is  state  vector  of  the  final  UCT  and  xaOb)  is  the  state  the  initial  UCT  state  propagated 
to  time  tg.  The  quiescent  propagated  state  is  calculated  by  taking  the  initial  condition  UCT 
and  propagating  it  under  the  assumption  that  no  maneuver  has  been  performed.  In  Euclidean 
distance,  only  the  positional  difference  is  considered  in  relating  the  two  states.  For  observation 
association  in  the  space  environment,  this  is  not  a  good  method  since  even  a  small  control  input 
can  cause  a  large  positional  diversion  over  time.  For  instance,  a  spacecraft  that  performs  a  plane 
change  maneuver  at  the  ascending  node  will  have  a  small  positional  difference  between  the 
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perturbed  and  homogenous  orbits  at  first,  but  after  a  quarter  of  an  orbit  period  the  out  of  plane 
displacement  will  be  significant.  Furthermore,  depending  on  the  time  of  observation,  objects  in 
very  different  orbits  reside  close  in  position.  In  each  of  these  cases,  a  geometric  distance  metric 
would  have  a  difficult  time  of  correctly  associating  uncorrelated  tracks. 

Alternately,  metrics  that  measure  the  distance  between  multivariate  normal  distributions, 
such  as  an  orbital  state  vector,  can  take  into  account  both  position  and  velocity  information  as 
well  as  uncertainty  in  the  state  vectors.  This  is  ideal  in  defining  a  metric  to  discriminate  between 
two  orbital  state  observations.  One  example  of  such  a  metric  is  the  Mahalonobis  distance,  which 
effectively  applies  a  weight  to  the  state  difference  along  each  axis  based  on  the  combined 
uncertainty  along  that  axis  in  the  observations  [10].  The  expression  for  the  Mahalonobis  distance 
is  given  in  Equation  10: 

dM(XA,PA.XB,PB,tB)  =  +  P  ,  (10) 

where  xb  and  Pb  are  the  full  state  and  covariance  information  of  the  state  estimate  from  the  final 
UCT,  andx^(fg^  and  Pa^b)  are  the  state  and  covariance  of  the  quiescent  propagated  state.  Similar 
to  Euclidean  distance,  the  propagated  homogenous  state  is  calculated  by  taking  the  initial 
condition  UCT  and  propagating  it  under  the  assumption  that  no  maneuver  has  been  performed, 
and  the  propagated  covariance  is  calculated  using  the  propagated  state  transition  matrix  to  update 
the  initial  covariance.  Mahalonobis  distance  improves  upon  Euclidean  distance  by  accounting  for 
uncertainty:  larger  uncertainty  in  any  axis  will  yield  a  lower  distance  metric  value.  A  value  of  1 
for  the  Mahalonobis  distance  is  a  1 -sigma  distance  (similarly,  a  value  of  2  is  a  2-sigma  distance, 
etc.).  Low  values  of  the  Mahalonobis  distance  yield  a  higher  confidence  in  correlation  of  the  full 
UCT,  accounting  for  both  the  full  state  (position  and  velocity)  and  state  uncertainty.  However, 
Mahalonobis  distance  is  still  at  its  core  a  state  differencing  metric.  Therefore,  it  still  has  the  issue 
in  orbital  mechanics  that,  depending  on  the  observation  times,  small  control  maneuvers  can 
cause  outsized  state  differences,  and  therefore  is  not  ideally  suited  for  orbital  mechanics 
applications.  Additionally,  the  Mahalanobis  distance  metric  does  not  characterize  maneuvers  -  it 
simply  accounts  for  differences  between  two  UCT  states. 

Both  the  Euclidean  distance  and  Mahalonobis  distance  metrics  are  calculated  alongside 
the  synthetic  and  WAAS  data  sensitivity  studies  (presented  in  section  4.2  and  4.3,  respectively). 
This  allows  a  direct  comparison  of  the  various  distance  metrics.  The  resulting  contour  plot  for 
Euclidean  distance  (in  kilometers)  is  shown  in  Figure  9.  First,  it  can  be  seen  that,  as  expected, 
increasing  SNR  does  not  affect  Euclidean  distance,  since  it  does  not  account  for  uncertainty.  The 
Euclidean  distance  increases  to  a  maximum  at  a  duration  of  one-half  orbit,  when  the 
homogenous  state  and  maneuvered  state  are  farthest  in  positional  separation.  The  metric  then 
decreases  back  toward  zero  at  an  orbit  duration,  when  the  orbits  cross  paths  at  the  line  of  nodes. 
This  behavior  is  cyclical  with  an  orbit  period.  The  difficulty  of  using  Euclidean  distance  for  UCT 
correlation  or  association  can  be  seen  as  the  value  of  the  metric  is  very  much  dependent  upon  the 
time  of  observation,  so  developing  thresholds  for  association  is  difficult.  For  instance,  the 
Euclidean  distance  breaks  down  as  an  association  metric  when  the  position  difference  is  small 
because  it  does  not  account  for  the  velocity  difference.  Even  RSOs  in  very  dissimilar  orbits  (such 
a  closed  orbit  and  a  fly-by  trajectory)  can  reside  in  close  proximity  for  a  short  time,  so  simply 
establishing  a  position  threshold  using  Euclidean  distance  is  not  sufficient. 

In  contrast,  the  Mahalonobis  distance  metric  plot,  shown  in  Figure  10,  is  much  more 
condensed  for  most  duration  values,  so  a  smaller  duration  time-frame  is  also  shown  in  Figure  10. 
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Here,  the  probability  volume  within  the  Mahalonobis  Distance  (MD)  isoprobability  surface  is 
computed  to  give  a  confidence  probability  that  the  final  state  is  consistent  with  the  initial  state, 
give  both  state  uncertainties.  It  should  be  noted  that  SNR  is  not  inherently  meaningful  for  use 
with  Mahalonobis  distance  since  minimum  detectable  cost  is  not  a  factor  in  the  equation  (the  ft 
matrix  from  Equation  4  ends  up  being  a  scaling  factor).  However,  for  better  comparison  between 
all  three  metrics,  Mahalonobis  distance  derived  probability  of  anomaly  was  also  plotted  against 
SNR.  At  durations  near  an  orbit  period,  the  Mahalonobis  distance  decreases,  showing  that  the 
UCTs  are  more  closely  correlated,  with  a  higher  likelihood  that  state  uncertainty  alone  explains 
the  observation  difference.  A  maneuver  is  harder  to  detect  at  these  points.  Conversely,  at  the  12- 
and  36-hour  durations,  where  the  out-of-plane  displacement  is  at  its  maximum,  the  Mahalonobis 
distance  is  considerably  higher,  showing  that  the  observed  state  change  likely  requires  a 
maneuver  or  other  anomaly.  Once  again,  this  is  due  to  the  inherent  state  differencing  in  the 
metric,  making  the  development  of  thresholds  difficult  for  application  to  arbitrary  UCT 
correlation. 


Figure  9.  Euclidean  Distance  Metric  (in  km)  for  Synthetic  Scenario  Sensitivity  Study 
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Figure  10.  Mahalonobis  Distance  Metric  for  Synthetic  Scenario  Sensitivity  Study 


The  advantage  that  the  control  distance  metric  has  over  these  other  distance  metrics  is 
that  it  is  able  to  account  for  the  specific  situation  of  orbital  mechanics  when  associating  UCTs. 
As  described  earlier,  even  small  control  perturbations  can  cause  a  large  state  divergence  over 
time,  even  though  the  orbits  are  actually  nearly  identical.  All  of  the  above  metrics  use  state 
differencing  to  compute  the  distance  between  the  observations,  whereas  the  control  distance 
metric  implemented  here  considers  the  optimal  control  required  to  connect  the  tracks.  In  the 
simplest  case,  an  un-maneuvered  object  should  have  a  control  distance  of  zero  between  any  two 
points  in  its  orbit,  providing  a  clear  indication  that  the  tracks  are  correlated.  In  this  case,  the 
geometric  or  normal  distribution  metrics  will  also  yield  zero  distance  metrics  because  the 
estimated  and  propagated  homogenous  states  are  identical.  For  maneuvered  tracks,  however,  the 
distance  metrics  described  in  this  section  do  not  provide  a  simple  convention  for  associating 
tracks,  as  it  would  need  to  involve  some  level  of  thresholding.  The  control  distance  metric 
described  in  this  study  provides  a  natural  method  for  correlating  maneuvered  tracks  by 
computing  the  maneuver  required  to  connect  the  tracks  and  accounting  for  uncertainty  to  test  the 
hypothesis  that  a  maneuver  occurred.  The  only  additional  logic  required  to  associate  tracks  is,  if 
the  probability  of  anomaly  is  high,  to  examine  the  history  of  the  maneuver  required.  If  the 
control  effort  required  exceeds  the  perceived  control  capability  of  the  RSO,  the  UCTs  can  be 
assumed  not  correlated,  and  vice  versa  for  a  reasonable  control  effort.  Therefore,  for  orbital 
mechanics  applications,  the  control  distance  metric  provides  an  intuitive  method  of  detecting 
maneuvers  and  correlating  UCTs. 

Other  distance  metrics,  such  as  Mahalonobis  distance,  are  commonly  used  in  UCT 
correlation,  and  perform  well  in  detecting  whether  two  UCTs  are  related  by  their  states  and 
uncertainties.  Mahalonobis  distance  can  be  useful  for  sensing  maneuvers  or  unmodeled 
dynamics,  as  the  distance  will  increase  drastically;  however,  it  cannot  provide  much  more 
information.  Therefore,  if  the  task  is  simply  to  detect  whether  a  maneuver  has  occurred, 
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Mahalonobis  distance  is  adequate.  In  comparison,  the  control  distance  metric  presented  here 
provides  an  anomaly  detection  capability  that  is  more  consistent  with  observation  duration  and 
takes  into  account  orbital  dynamics.  Additionally,  since  the  nominal  trajectory  must  be 
computed,  the  trajectory  connecting  the  UCTs  is  also  obtained  in  the  process,  which  can  be 
reviewed  to  inspect  realism  and  feasibility  of  the  anomaly  hypothesis. 

The  differences  between  the  metrics  lend  themselves  to  different  applications.  Consider  a 
scenario  where  three  space  objects  are  detected  in  close  proximity.  At  a  later  observation,  three 
objects  are  also  detected,  but  some  or  all  of  them  have  modified  orbit  element  slots.  Mahalonobis 
distance  could  be  used  to  hypothesize  whether  the  objects  had  indeed  maneuvered,  based  on  state 
differences  and  uncertainty,  but  it  merely  provides  a  statistical  similarity  set  of  the  data.  For 
instance,  Mahalonobis  distance  cannot  account  for  maneuvering  objects.  With  control  distance 
metrics,  the  same  information  and  more  can  be  gathered.  The  trajectories  required  to  connect 
each  UCT  can  be  calculated,  providing  an  ordered  pair  of  the  most  likely  correlations  and  their 
respective  trajectory  or  maneuver  requirements.  In  the  three-RSO  scenario,  the  UCT 
combinations  that  yield  the  lowest  total  control  cost  would  be  the  most  likely  correlations,  and 
the  connecting  trajectories  can  be  reviewed  to  confirm  the  hypotheses. 

In  addition  to  comparisons  between  metrics  with  synthetic  data,  the  metrics  were  used  to 
investigate  maneuver  detection  in  the  WAAS  data.  The  Euclidean  distance  metric  plot,  shown  in 
Figure  11,  bears  resemblance  to  the  synthetic  data  plot;  the  oscillations  are  smaller  in  magnitude 
in  this  case,  but  the  Euclidean  distance  does  show  periodic  behavior  with  the  orbit  period.  The 
differences  between  Figure  9  and  11  can  be  attributed  to  the  unmodeled  dynamics  and  the 
spacecraft’s  continuous  maneuvering  to  reject  those  perturbations. 

The  Mahalonobis  distance  plot,  in  Figure  12,  also  bears  resemblance  to  the  synthetic  data 
plot  in  Figure  10,  showing  the  same  improved  anomaly  detection  capability  (higher  Mahalonobis 
distance  for  a  given  SNR)  at  the  12-  and  36-hour  durations,  where  the  out-of-plane  displacement 
difference  is  greatest.  The  generally  lower  values  of  Mahalonobis  distance  at  other  durations  are 
also  likely  due  to  unmodeled  dynamics  and  the  continuously  thrusting  behavior  of  the  spacecraft 
data.  Mahalonobis  distance  is  only  able  to  detect  that  a  maneuver  has  occurred,  but  control 
distance  can  determine  more  information;  for  instance,  the  type  of  maneuver  can  be  determined 
based  on  the  control  and  trajectory  from  the  trajectory  optimization  step.  In  contrast,  control 
distance  can  both  confirm  maneuvers  and  provide  the  most  likely  combination  of  states, 
providing  added  data  over  the  existing  metrics. 

On  the  surface,  it  may  appear  that  the  consistently  lower  SNR  (or  equivalently,  track  state 
uncertainty)  required  to  associate  UCTs  would  lead  a  prudent  SSA  operator  to  use  the 
Mahalanobis  distance  metric.  However,  based  on  the  control  distance  metric  results,  it  is  the 
view  of  the  authors  that  such  reliance  on  MD  data  association  may,  at  particular  points  in  the 
orbit  or  for  specific  durations,  prematurely  associate  data  and  eliminate  potentially  useful 
association  hypotheses.  Said  differently,  the  Mahalanobis  distance  metric  may  lead  to  missed 
associations.  This  is  true  namely  because  the  Mahalanobis  distance  metric  assumes  that  no 
maneuvers  have  taken  place.  The  control  metric  distribution  given  the  same  information 
necessarily  requires  lower  uncertainty  (higher  SNR)  to  claim  that  an  anomaly  has 
occurred  because  the  assumed  system  may  include  controlled  thrust. 
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Figure  11.  Euclidean  Distance  Metric  (in  km)  for  WAAS  Scenario  Sensitivity  Study 


Figure  12.  Mahalonobis  Distance  Metric  for  WAAS  Scenario  Sensitivity  Study 


An  observation  of  the  Mahalanobis  distance  computation  method  must  be  made.  As 
previously  stated,  the  MD  implementation  for  this  report  does  not  model  maneuvers  or 
perturbations.  One  way  to  model  such  perturbations  is  to  include  process  noise  (dynamics 
uncertainty).  Then,  as  the  uncertainty  covariance  is  propagated,  the  volume  of  the  uncertainty 
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necessarily  increases,  and  Mahalanobis  distances  for  all  associations  naturally  reduce.  Thus, 
when  process  noise  is  included,  fewer  track  association  hypotheses  are  eliminated,  as  in  the 
control  distance  results.  The  problem  with  this  approach,  however,  is  two-fold:  1)  maneuvers 
result  in  structured  accelerations,  and  are  not  well  modeled  with  random  accelerations  assumed 
using  process  noise,  and  2)  to  appropriately  choose  the  process  noise  to  capture  the  maneuvering 
trajectory,  one  must  first  know  the  approximate  size  of  the  maneuver,  resulting  in  a  ‘chicken  or 
the  egg’  problem. 

To  summarize,  the  Mahalanobis  distance  metric  functions  well  as  an  UCT  association 
method  when  space  objects  (SOs)  are  quiescent  or  when  there  are  no  other  nearby  SOs. 
However,  if  maneuvers  occur,  particularly  when  clusters  of  SOs  are  concerned,  the  Mahalanobis 
distance  metric  eliminates  credible  association  hypotheses  because  active  structured  maneuvers 
are  not  considered. 

5.0  CONCLUSIONS 

This  research  demonstrates  the  application,  effectiveness,  and  limitations  of  the 
DAMMM  algorithm.  The  algorithm  is  designed  to  solve  the  two-point  boundary  value  problem 
of  space  object  observations  to  address  the  needs  of  object  track  correlation  and  maneuver 
characterization.  The  method  makes  use  of  control  cost  as  a  distance  metric,  constructing  a 
trajectory  to  connect  the  boundary  conditions  under  Keplerian  dynamics  and  optimal  control 
assumptions.  A  scenario  is  constructed  to  mirror  a  station-keeping  maneuver,  and  the  results 
were  compared  to  actual  station-keeping  maneuvers  from  WAAS  observational  data  of  the 
Galaxy  1 5  spacecraft.  Additionally,  the  limitations  of  the  algorithm  are  tested  to  determine  the 
sensitivity  of  the  algorithm  to  sensor  noise  and  observation  times.  In  general,  anomaly 
probability  is  consistent  with  varying  duration  and  constant  SNR.  Finally,  the  control  distance 
metric  is  compared  to  other  common  distance  metrics  to  highlight  the  advantages  of  using 
control  distance  to  correlate  object  tracks. 
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APPENDIX  A:  DAMMM  Software  Framework  Code  Description 


This  appendix  describes  the  DAMMM  software  framework  as  shown  in  Figure  1  of  the 
report.  Each  major  component  of  the  framework  is  described  with  the  inputs  required  and 
outputs  provided.  A  script  can  be  written  as  a  wrapper  to  handle  the  input  tracks  and  tie  the 
functions  together  as  shown  in  Figure  1,  tailored  to  the  specific  scenario  encountered.  All  code  is 
implemented  in  MATLAB. 


Trajectory  Optimization  -  trajOptimizer.m 


[optTraj]  =  traj Optimizer ( scenario, settings ) 

Solves  a  two-point  boundary  value  problem  to  corresponding  to  the  input 
scenario.  Using  Keplerian  orbital  dynamics,  the  optimizer  seeks  to  find 
the  thrust  accelerations  required  to  optimally  maneuver  a  spacecraft 
to  match  the  boundary  conditions  using  a  quadratic  cost  function, 
consistent  with  variable-specific  impulse  engine  performance.  The 
constrained  minimization  function  (fmincon)  utilizes  well-known 
gradients  for  the  cost  function  and  dynamics  constraints  to  improve 
convergence  performance,  and  the  settings  structure  input  can  be  used  to 
further  tune  optimizer  performance. 


sINPUT 
5  scenario 


settings 


structure  defining  scenario,  must  contain  the  following 


,  km/ s ) 
km/ s ) 


structure  containing  settings  for  optimization  routine, 
contains  the  following  fields: 

' maxFuncCalls '  -  1x1  scalar,  maximum  number  of  steps  to 
use  in  optimization  routine 


o 

o 

fields : 

o 

o 

'xO  ' 

-  6x1 

vector. 

initial  ECI  state  (km 

o 

o 

'xf ' 

-  6x1 

vector. 

final  ECI  state  (km. 

o 

o 

'to ' 

-  1x1 

scalar. 

initial  time  (s) 

o 

o 

'  tf ' 

-  1x1 

scalar. 

final  time  (s) 

o 

o 

'dt ' 

-  1x1 

scalar. 

time  step  (s) 

default:  3000 

1x1  scalar,  absolute  convergence 
tolerance  in  optimizer 
1x1  scalar,  relative  convergence 
tolerance  in  optimizer 

1x1  scalar,  termination  tolerance  for 
changes  in  decision  variable 
default:  le-10 

1x1  scalar,  tolerance  on  constraint 
violation 
default:  le-6 

boolean,  true  to  generate  a  guess  of 
control  inputs  to  prime  the  optimizer 
' useGradients '  -  boolean,  true  to  use  gradients  in  cost 
function  and  trajectory  constraints, 
default:  'true' 

string,  selects  display  mode  for 
optimizer,  review  fmincon  'Display' 
for  options 


'  abstol ' 


' reltol ' 


' tolx ' 


' tolcon ' 


' genGuess ' 


' displayMode ' 
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%OUTPUT 
%  optTraj 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o, 

o 


structure  containing  resulting  optimal  trajectory  data 


' time ' 

-  Ixn  vector,  simulation  time  steps  (s) 

' tra j  ectory ' 

-  6xn  vector,  state  at  each  time  step  for 
optimal  trajectory  (km,  km/s) 

' control ' 

-  3xn  vector,  thrust  accelerations  at  each 
time  step  for  optimal  trajectory 
(km/s"2) 

'dV 

-  1x1  scalar,  delta-v  computed  by 

integrating  thrust  acceleration  output 
( km/ s ) 

' cost ' 

-  1x1  scalar,  optimal  quadratic  control 
cost  (km'^2/s^3) 

Fit  Adjoints  -  fitAdjoints.m 


o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 


[xstar , pstar , ustar , tstar ]  =  f itAdj oints (xinp, ulnp, tinp, mu, Re) 

Convert  the  direct  trajectory  (x(t),u(t))  with  quadratic  control  cost  to 
an  equivalent  indirect  trajectory  (e.g.,  xO,  pO) .  Uses  input  control 
array  to  generate  costate  guess  under  assumptions  of  optimal  control, 
iteratively  update  costate  guesses  until  costate  guess  converges. 
Propagates  new  indirect  trajectory  using  Keplerian  dynamics  to  check 
errir  between  new  optimal  trajectory  and  input  trajectory. 


INPUT 

xinp  -  6xn  array,  direct  trajectory  states  (km,  km/s) 
ulnp  -  3xn  array,  direct  trajectory  controls  (km/s"^2) 
tinp  -  Ixn  vector,  direct  trajectory  times  (s) 
mu  -  1x1  scalar,  gravitational  parameter  of  planet 
Re  -  1x1  scalar,  equatorial  radius  of  planet 


OUTPUT 

xstar  -  6xn  array,  optimal  indirect  trajectory  states  (km,  km/s) 

pstar  -  6xn  array,  optimal  indirect  trajectory  costates 

ustar  -  3xn  array,  optimal  indirect  trajectory  controls  (km/s'^P) 

tstar  -  Ixn  vector,  optimal  indirect  trajectory  times  (s) 


Compute  Cost  Distribution 

g _ 

o - 

%  [pdf_cost ,  quad_cost ,  min_cost ,  mu_Pu,  sig_Pu] 

%  =  genCostDist (xvec, pvec, tvec, BCs , Pkk, params ) 

o 

o 

%  Generate  cost  distribution  for  nominal  indirect  trajectory  given  boundary 
%  condition  and  systemic  uncertainty  inputs . 

o 

o 

%INPUTS 

%  xvec  6xn  km,  km/s,  mean  state  space  trajectory  over  t 

%  pvec  6xn  (),  ()/s,  mean  adjoint  trajectory  over  t 
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o\o  o\o 


tvec  Ixn  s,  time  vector 

BCs  Cell  array  of  boundary  conditions,  i=l,...,Np.  Each  BC  has  four 

parameters 

%  BC{i}.t  1x1  scalar  time  of  boundary  condition 

%  BC{i}.P  pxp  covariance  matrix  of  ith  boundary  condition  uncertainty 
%  BC{i}.y  pxl  nominal  (mean)  ith  boundary  condition  value 

%  BC{i}.H  pxn  partial  derivative  of  ith  boundary  condition  (\partial  g 

/partial  x) 

%  Pkk  2x2  covariance  matrix  of  dynamics  parameter  uncertainty ( ies ) 

%  (atmospheric  density  and  ballistic  coefficient) 

%  params  Cell  array  of  dynamics  parameters  parameters 


o 

o 

params { 1 } 

1x1 , 

km^2/s^3,  MU,  Central  body  gravitational  constant 

o 

o 

params { 2 } 

1x1 , 

0,  J2  coefficient 

o 

o 

params { 3 } 

1x1 , 

0,  J3  coefficient 

o 

o 

params { 4 } 

1x1 , 

km,  Radius  of  central  body 

o 

o 

params { 5 } 

1x1 , 

kg/m^2.  Atmospheric  density 

o 

o 

params { 6 } 

1x1 , 

kg/m^2.  Ballistic  coefficient 

o 

o 

params { 7 } 

1x1 , 

km/s^2.  Maximum  available  thrust  acceleration 

o 

o 


%OUTPUTS 

%  pdf_cost  -  Ixr  probability  distribution  of  quadratic  control  costs 

%  cost_dist  -  Ixr  quadratic  cost  values  for  pdf_cost  (km^2/s^3) 

%  min_cost  -  1x1  minimum  detectible  quadratic  cost  (expectation  of 
%  systemic  uncertainty)  (km'^2/s^3) 

%  mu_Pu  -  1x1  first  moment  of  cost  distribution 

%  sig_Pu  -  1x1  second  moment  of  cost  distribution 

o 

o 

o, 

o 


Test  Anomaly  Hypothesis  -  det_anomaly _cost.m 


g _ 

o - 

% [p_anomaly, details ]  =  det_anomaly (pdf_DV, DV, method, det_params ) 

o 

o 

%  Determine  the  probability  that  an  anomaly  has  occurred.  Given  input 
%  distribution  and  trajectory  data,  determine  probability  that  the  observed 
%  state  change  is  attributable  to  anomaly  or  uncertainty  in  state. 

o 

o 

%INPUT 

%  pdf_cost  -  Ixr  probability  distribution  of  Delta  V 
%  cost_dist  -  Ixr  Delta  V  values  for  pdf_DV 

%  method  string  argument  for  detection  method  to  be  used.  Options 

%  are  'HSA'  (Holzinger-Scheeres-Alf riend) 

%  det_params  -  cell  of  detection  parameters  for  the  method  selected 
%  mu_Pu  -  1x1  first  moment  of  cost  distribution 

%  sig_Pu  -  1x1  second  moment  of  cost  distribution 

o 

o 

%OUTPUT 

%  p_anomaly  -  1x1  probability  that  an  anomaly  has  occured 
%  details  -  method  dependent  structure  of  ancillary  information 

o 

o 

o. 

o 
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APPENDIX  B:  Theory  Extension 

This  appendix  describes  novel  theoretical  results  that  may  be  used  to  compute  the 
distribution  of  quadratic  control  distances  using  Gaussian  mixtures  as  boundary  condition 
uncertainties. 

The  quadratic  control  cost  for  a  single  term  of  the  Gaussian  sum  connecting  Gaussian 
initial  distribution  /  to  Gaussian  final  distribution  j  may  be  described  as 

P,j  =  kj  +  mJjVij  +  VijNijVij  (B-1) 

where  lij  E  R,  mfj  E  Nij  E  and  the  random  variable  V  ~  N{jiij,Py  ij),  where 

Defining  a  new  random  variable  Wij 

such  that 

where  now  W  ~  Al(0,  Py.tj)-  The  cost  Pij  can  be  re-written  in  terms  of  this  new  random  variable 
using  direct  substitution: 

^ij  ~  hj  T  T  T  {_Pij  T  ^t;)  T 

=  hj  +  +  mJjWij  +  pJjNijfiij  +  +  WjjNijWij 

=  {lij  +  +  iifjN  ijiiij)  +  (niij  +  2^iJjNij)Wij  +  WjjNijWij 


If  the  following  definitions  are  made: 


hj  =  hj+  +  pJjNijUij 

^ij  =  mjj  +  2fiJjNij 


then  Equation  B-1  may  be  rewritten  in  a  nearly  identical  form  using  the  zero-mean  random 
variable  Wij  as 


Pt,  =  kj  +  +  WJiNtjWt,  (B-2) 

From  Holzinger  et  al.  (see  reference  8),  Appendix  B,  the  analytic  first  and  second  moments  of 
Pij  may  be  written  as: 

E[Pjy]  —  Pp^ij  hj  "f  ^ijPij  "f  Pij^ijPij  "f  Tr[lViyPy 

~  ^P,ij  ~  iv^ij  T  2PiylVjy)  Py  ij{Ttlij  -f  2pijNij^  +  2Tv[NijPy  ijNijPy  ij^ 

The  total  cost  of  all  i  initial  boundary  conditions  and  j  final  boundary  conditions  is  the  weighted 
sum  of  the  individual  costs  between  each  i  and  j: 
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i  j 

where  W;  is  the  weight  of  the  f’'  initial  boundary  condition  and  Wj  is  the  weight  of  the  /*  final 
boundary  condition.  Thus,  the  analytic  expected  value  of  the  initial  and  final  Gaussian  sum 
boundary  conditions  may  be  written  as: 

E[P]  =  E 

^  j 

=  '^^WiWjE[Pij] 

i  j 

i  j 

which  leads  to  the  final  result: 

Up  =  ZiZ;  +  rnjjfiij  +  fi-jNijfiij  +  Tr[NijPv,ij])  (B-3) 
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LIST  OF  SYMBOLS,  ABBREVIATIONS,  AND  ACRONYMS 


AFRL 

CDF 

DAMMM 

ECEF 

GEO 

MD 

PDF 

PI 

RIC 

RAM 

RSO 

SNR 

SO 

SSA 

Sow 

TPBVP 

UCT 

VSI 

WAAS 


Air  Force  Research  Laboratory 
Cumulative  distribution  function 

Data  association  using  a  minimum-fuel  maneuver  metric 

Earth-centered  Earth-fixed 

Geostationary  Earth  orbit 

Mahalonobis  Distance 

Probability  density  function 

Principal  investigator 

Radial,  in-track,  and  cross-track 

Random  access  memory 

Resident  space  object 

Signal-to-noise  ratio 

Space  object 

Space  situational  awareness 
Statement  of  work 

Two-point  Boundary  Value  Problem 
Uncorrelated  track 
Variable  specific  impulse  engines 
Wide  Area  Augmentation  System 
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